Within this chapter developmental plasticity was explored within Acanthochromis polyacanthus that were collected from two different regions (i.e., low-latitude, Cairns, and high-latitude, Mackay). Fish were held in common garden experiments at 28.5 C. Clutch data was collected along with parental morphometeric data to determine how fish from each population performed at 28.5 C, a temperature that was shown to produce similar absoluate aerobic scope performances in a previous study [link]. Once hatched offspring were placed into three different temperature treatments, 28.5, 30, and 31.5 C. After approximately 5-6 months offspring length and weight was measured, as well as CTmax and respiration during CTmax trials.
library(tidyverse) # data manipulation
library(ggpubr) # figure arrangement
library(brms) # Bayesian models
library(StanHeaders)# needed to run Bayesian models
library(rstan) # needed to run Bayesian models
library(standist) # needs to be installed
library(bayesplot) # needed for MCMC diagnostics
library(DHARMa) # model validation
library(ggdist) # partial plots
library(tidybayes) # partial plots
library(broom.mixed) # model investigation
library(emmeans) # pairwise comparisons
library(rstanarm) # pairwise comparisons - need for emmeans
library(ggpubr) # arranging figures
library(mgcv) # GAM
library(modelr) # plottingThis data set has no point outliers however, in the next chunk 3 samples will be discarded due to poor data quailty
rows_of_data <- count(lat_resp_dat, sampleID)
lat_resp_dat2 <- lat_resp_dat |>
mutate(dev.temp = as.factor(dev.temp),
replicate = str_sub(sampleID, -1,-1),
population = factor(population)) |>
# number of observations = 5758
filter(sampleID != "56.CARL.137.28,5.1", # 5777 - 76 = 5682
sampleID != "56.CARL.137.28,5.2", # 5701 - 64 = 5618
sampleID != "60.LCKM.152.30.1" # 5637 - 76 = 5542
) |>
filter(time_lag_sec >2001) |> # remove samples from first 5 cycles (i.e., first 33 minutes)
group_by(sampleID) |>
mutate(max_value_index = which.max(rate.output),
row_number = row_number(),
max.rate = max(rate.output),
low.rate = mean(rate.output[order(rate.output)[1:10]]),
net.rate = max.rate - low.rate) |>
filter(row_number <= max_value_index) |>
ungroup() |>
mutate(region = factor(region),
dev.temp =factor(dev.temp),
level = as.factor(case_when(tank >= 1 & tank <= 199 ~ 1,
tank >= 200 & tank <= 299 ~ 2,
tank >= 300 & tank <= 399 ~ 3,
TRUE ~ NA_real_)),
female = factor(female))eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) +
geom_point(alpha=0.5, color = "black") +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
theme_classic() +
ggtitle("All data points - 2nd order polynomial")
eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ dev.temp) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
nrow = 6,
ncol=1); eda.figeda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) +
geom_point(alpha=0.5, color = "black") +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE), color = "red") +
theme_classic() +
ggtitle("All data points - 2nd order polynomial")
eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE), color = "red") +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE)) +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE)) +
facet_wrap(~ dev.temp) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE), color = "red") +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE)) +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
nrow = 6,
ncol=1); eda.figFirst we will place random factors only within out model and see if they do a good job explaining the variance within our model. We will also be include priors which will be based off of out length data.
Hypothesis test will include:
lat_resp_dat2 |>
group_by(region,dev.temp) |>
summarise(mean(na.omit(rate.output)),
sd(na.omit(rate.output))) ## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.
f.model.null <- bf(rate.output ~ 1,
family = gaussian())
model.null <- brm(f.model.null,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
f.model1 <- bf(rate.output ~ 1 + (1| female) + (1| tank),
family=gaussian())
model1 <- brm(f.model1,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
f.model2 <- bf(rate.output ~ 1 + (1| female) + (1| tank) + (1| level),
family=gaussian())
model2 <- brm(f.model2,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## Warning: There were 13 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
f.model3 <- bf(rate.output ~ 1 + (1| female) + (1| tank) + (1| level) + (1| population),
family=gaussian())
model3 <- brm(f.model3,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## Warning: There were 5 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
f.model4 <- bf(rate.output ~ 1 + (1| female) + (1 | tank) + (1| level) + (1| population)+ (1| clutch_order),
family=gaussian())
model4 <- brm(f.model4,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
## Output of model 'model.null':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -31469.1 78.1
## p_loo 2.5 0.4
## looic 62938.2 156.2
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.0]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model1':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30662.6 99.6
## p_loo 57.6 3.1
## looic 61325.1 199.1
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model2':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30662.7 99.5
## p_loo 57.6 3.0
## looic 61325.4 198.9
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model3':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30662.7 99.6
## p_loo 58.1 3.2
## looic 61325.5 199.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model4':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30662.5 99.4
## p_loo 57.5 3.0
## looic 61325.0 198.8
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## model4 0.0 0.0
## model1 -0.1 0.4
## model2 -0.2 0.4
## model3 -0.2 0.4
## model.null -806.6 48.2
lat_resp_dat2 |>
group_by(region,dev.temp) |>
summarise(mean(na.omit(rate.output)),
sd(na.omit(rate.output))) ## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.
rate.priors <- prior(normal(430, 0.25), class="Intercept") +
prior(normal(0, 40), class="b")
prior(student_t(3, 0, 195), class = "sigma")f.model1.p1 <- bf(rate.output ~ 1 +
dev.temp*region + Temperature +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.p1 <- brm(f.model1.p1,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
#--- distribution check ---#
pp_check(model1.p1, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.p1, ndraws=250, summary=FALSE)
model1.p1_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.p1_resids) ; testDispersion(model1.p1_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98596, p-value = 0.704
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
f.model1.p2 <- bf(rate.output ~ 1 +
dev.temp*region + poly(Temperature, 2) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.p2 <- brm(f.model1.p2,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
#--- distribution check ---#
pp_check(model1.p2, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.p2, ndraws=250, summary=FALSE)
model1.p2_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.p2_resids) ; testDispersion(model1.p2_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98971, p-value = 0.688
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
f.model1.p3 <- bf(rate.output ~ 1 +
dev.temp*region + poly(Temperature, 3) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.p3 <- brm(f.model1.p3,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
#--- distribution check ---#
pp_check(model1.p3, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.p3, ndraws=250, summary=FALSE)
model1.p3_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.p3_resids) ; testDispersion(model1.p3_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98512, p-value = 0.616
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
## Output of model 'model1.p1':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -29946.0 94.1
## p_loo 60.3 3.1
## looic 59892.1 188.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model1.p2':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30525.1 100.7
## p_loo 58.9 3.2
## looic 61050.1 201.3
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.2]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model1.p3':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30516.1 100.1
## p_loo 58.7 3.1
## looic 61032.1 200.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.2]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## model1.p1 0.0 0.0
## model1.p3 -570.0 31.3
## model1.p2 -579.0 31.5
f.model1.gamk4 <- bf(rate.output ~ 1 +
t2(dev.temp,region,Temperature, k=4, bs="fs") +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank),
family=gaussian())
model1.gamk4 <- brm(f.model1.gamk4,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
#--- distribution check ---#
pp_check(model1.gamk4, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.gamk4, ndraws=250, summary=FALSE)
model1.gamk4_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.gamk4_resids) ; testDispersion(model1.gamk4_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.97732, p-value = 0.464
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
### {-}
f.model1.gamk3 <- bf(rate.output ~ 1 +
t2(dev.temp,region,Temperature, k=3, bs="fs") +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank),
family=gaussian())
model1.gamk3 <- brm(f.model1.gamk3,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
#--- distribution check ---#
pp_check(model1.gamk3, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.gamk3, ndraws=250, summary=FALSE)
model1.gamk3_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.gamk3_resids) ; testDispersion(model1.gamk3_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9811, p-value = 0.536
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
### {-}
## Warning in AIC.default(model1.gamk3, model1.gamk4): models are not all fitted
## to the same number of observations
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: rate.output ~ 1 + t2(dev.temp, region, Temperature, k = 4, bs = "fs") + scale(mass, center = TRUE, scale = TRUE) + (1 | female) + (1 | tank)
## Data: lat_resp_dat2 (Number of observations: 4776)
## Draws: 2 chains, each with iter = 5000; warmup = 500; thin = 5;
## total post-warmup draws = 1800
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sds(t2dev.tempregionTemperature_1) 809.91 168.24 551.16 1213.34 1.00
## sds(t2dev.tempregionTemperature_2) 878.19 186.51 593.72 1318.30 1.00
## Bulk_ESS Tail_ESS
## sds(t2dev.tempregionTemperature_1) 1300 1719
## sds(t2dev.tempregionTemperature_2) 1681 1669
##
## Group-Level Effects:
## ~female (Number of levels: 15)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 113.90 32.16 62.97 185.66 1.00 1535 1705
##
## ~tank (Number of levels: 52)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 82.15 10.73 64.41 107.17 1.00 1721 1643
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 430.00 0.26 429.48 430.50 1.00
## scalemasscenterEQTRUEscaleEQTRUE 19.84 3.65 12.56 27.16 1.00
## Bulk_ESS Tail_ESS
## Intercept 1766 1635
## scalemasscenterEQTRUEscaleEQTRUE 1894 1643
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 114.30 1.18 111.91 116.61 1.00 1713 1567
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#model1.re.wo |> get_variables()
model1.gamk4 |> gather_draws(`b_.*|sigma`, regex =TRUE) |>
median_hdci()model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="dev.temp") |> summary(infer=TRUE)model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="region") |> confint()## region dev.temp emmean lower.HPD upper.HPD
## core 28 407 325 502
## leading 28 404 297 516
## core 30 447 367 520
## leading 30 387 268 494
## core 31 452 360 531
## leading 31 468 374 580
##
## Point estimate displayed: median
## HPD interval probability: 0.95
model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="region") |> confint()## $emtrends
## dev.temp region Temperature.trend lower.HPD upper.HPD
## 28 core 42.04 29.76 56.2
## 30 core 1.39 -10.50 12.5
## 31 core 20.20 8.77 31.8
## 28 leading -1.66 -15.93 11.0
## 30 leading 1.03 -13.20 14.5
## 31 leading 4.25 -8.60 16.4
##
## Point estimate displayed: median
## HPD interval probability: 0.95
##
## $contrasts
## contrast estimate lower.HPD upper.HPD
## dev.temp28 core - dev.temp30 core 40.607 24.078 60.12
## dev.temp28 core - dev.temp31 core 21.926 3.890 39.59
## dev.temp28 core - dev.temp28 leading 44.028 24.126 61.31
## dev.temp28 core - dev.temp30 leading 41.138 21.869 60.16
## dev.temp28 core - dev.temp31 leading 37.628 20.663 57.08
## dev.temp30 core - dev.temp31 core -18.947 -35.413 -2.84
## dev.temp30 core - dev.temp28 leading 2.996 -14.071 20.71
## dev.temp30 core - dev.temp30 leading 0.223 -17.503 19.71
## dev.temp30 core - dev.temp31 leading -3.028 -19.615 12.87
## dev.temp31 core - dev.temp28 leading 21.958 4.079 38.88
## dev.temp31 core - dev.temp30 leading 19.222 -0.815 36.84
## dev.temp31 core - dev.temp31 leading 15.969 -1.138 32.29
## dev.temp28 leading - dev.temp30 leading -2.709 -23.119 15.26
## dev.temp28 leading - dev.temp31 leading -5.889 -23.118 13.74
## dev.temp30 leading - dev.temp31 leading -3.028 -20.894 16.54
##
## Point estimate displayed: median
## HPD interval probability: 0.95
resp.plot <- lat_resp_dat2 |>
group_by(region, dev.temp) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_fitted_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=Temperature, color=exp_group)) +
geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) +
theme_classic() +
facet_wrap(~dev.temp); resp.plotresp.plot2 <- lat_resp_dat2 |>
group_by(region, dev.temp) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_fitted_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=Temperature, color=exp_group)) +
geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) +
theme_classic() +
facet_wrap(~region); resp.plot2resp.plot3 <- lat_resp_dat2 |>
group_by(region, dev.temp) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_fitted_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=Temperature, color=exp_group)) +
geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) +
theme_classic(); resp.plot3